Modeling the Simultaneous Effects of Particle Size and Porosity in Simulating Geo-Materials

The particle discrete element method (PDEM) is widely used to simulate rock and soil materials to obtain stress and strain. However, there are three shortcomings: (1) Single sphere or ellipsoids directly replace the soil particles; (2) it treats the diameters of spheres or ellipsoids as the soil particle size; (3) the overlapping particle volume is not deducted in calculating the porosity. Hence, it is difficult for the simulation of the geological body to agree with reality. This research found a rotation calculation model and a pixel counting method to make joint soil particles more accurately simulate geological materials to solve the three shortcomings. The model successfully obtained the gradation curve and porosity of the simulated geological body with joint particles. This research will further enrich and broaden the application prospects of PDEM and provide a reference for scientific research and engineering fields in geological engineering, geotechnical engineering, and petroleum engineering.


Introduction
In order to obtain the mechanical properties of rock and soil mass, a large number of in situ and laboratory tests are carried out during the construction of hydraulic and geotechnical engineering. The cost of many in situ and laboratory tests is high. After the onsite sampling, the indoor test cannot guarantee the sample's original stress; the test results largely deviate from the actual one. Simultaneously, different lithologies and sampling locations make it impossible to apply these test results widely. The numerical calculation method can yield geomechanical parameters, including stress-strain and porosity at specific depths under these standard preconditions to achieve these goals: low costs, maintaining the original stress and simulating different soils and rock. Some methods can model geological materials, mainly the finite element method (FEM) and the discrete element method (DEM). The finite element method has solved many theoretical, technical, and engineering problems in geology [1,2]. However, many issues still cannot be solved. First, the finite element method is unsuitable for large deformation [3]. Second, the discontinuous rock mass cannot use the technique, as these bodies contain voids, cracks, and fissures. Third, this method cannot directly obtain the rock and soil porosity. The following formula containing stress-strain parameters calculates the porosity. The strain and stress result from the finite element method of analyzing soils/rock. Then, the new porosity is Equation (1). where n new = e 0 −∆ 1+e 0 −∆ , n 0 = e 0 1+e 0 , θ t = 1−2υ E Θ, θ t = ∆ 1+e 0 , n 0 is the original porosity, E is the elastic modulus, υ is the Poisson's ratio, e 0 is the original void ratio, ∆ is the pore change/reduction, e 0 is the original void ratio, Θ = σ x + σ y + σ z , σ x , σ y , and σ z are the stress in the x, y, and z direction, θ t is the volumetric strain, and n new is the new porosity. The void is pressed and reduced, but the volume of the soil/rock particles does not diminish. So, the new porosity is calculated by the stress/strain. The method is not direct.
The discrete element method was proposed by Cundall and Strack [4] and was widely used in analyzing large rock and soil deformation. The method can solve significant deformation problems simulate a discontinuous body and granular soils [4]. Furthermore, soils contain many cracks, faults, micropores, and microchannels and are not continuous or uniform bodies. There are many mechanical interactions between cracks and particles [5,6]. The discrete element method is widely used in geotechnical, geological, and hydraulic fields [7,8].
Specifically, the discrete element method has some branches. For example, the commercial software of the sphere discrete element method [9][10][11][12] is PFC and Yade, esys-particle and LIGGGHTS ( Figure 1); the commercial software of the discontinuous medium (jointed rock mass) discrete element method [13][14][15] is UDEC/3DEC, FLAC and EDEM, the discrete lattice points discrete element method [2,16], such as Xsite. However, the particle discrete element method (PDEM) is more widely used in materials processing. , n0 is the original porosity, E is the elastic modulus, υ is the Poisson's ratio, e0 is the original void ratio,  is the pore change/reduction, e0 is the original void ratio, Θ = σx + σy + σz, σx, σy, and σz are the stress in the x, y, and z direction, θt is the volumetric strain, and nnew is the new porosity. The void is pressed and reduced, but the volume of the soil/rock particles does not diminish. So, the new porosity is calculated by the stress/strain. The method is not direct. The discrete element method was proposed by Cundall and Strack [4] and was widely used in analyzing large rock and soil deformation. The method can solve significant deformation problems simulate a discontinuous body and granular soils [4]. Furthermore, soils contain many cracks, faults, micropores, and microchannels and are not continuous or uniform bodies. There are many mechanical interactions between cracks and particles [5,6]. The discrete element method is widely used in geotechnical, geological, and hydraulic fields [7,8].
Specifically, the discrete element method has some branches. For example, the commercial software of the sphere discrete element method [9][10][11][12] is PFC and Yade, esys-particle and LIGGGHTS ( Figure 1); the commercial software of the discontinuous medium (jointed rock mass) discrete element method [13][14][15] is UDEC/3DEC, FLAC and EDEM, the discrete lattice points discrete element method [2,16], such as Xsite. However, the particle discrete element method (PDEM) is more widely used in materials processing. Polarized microscopic marble slice and crystal particle simulation (a) Single-particle/circular particle/ball particle; (b) one crystal that is a joint particle/crystal particle; (c) partial enlargement of a joint particle; (d) crystal in marble slice. 1, 2, 3, 4, and 5 are five crystals in a marble slice, this scale is approximate scale. The research on the particle discrete element method (PDEM) [2,9,17] has achieved fruitful results, such as the reference parameters of the Brazilian disc failure test and PFC2d simulation [9], the laws of various mechanical parameters, the coupling relationship with fluids [12,18], the relationship between the particle size and equilibrium position for the particles in the fluid channel [19], the deformation characteristics of geological bodies, and disaster simulations [20], there are many mechanical interactions between cracks and the particles [5,6]. The method measures particle velocity in water and multiphase flow by particle image velocimetry [21,22]. These studies simulate the soil's circular Figure 1. Polarized microscopic marble slice and crystal particle simulation (a) Single-particle/ circular particle/ball particle; (b) one crystal that is a joint particle/crystal particle; (c) partial enlargement of a joint particle; (d) crystal in marble slice. 1, 2, 3, 4, and 5 are five crystals in a marble slice, this scale is approximate scale. The research on the particle discrete element method (PDEM) [2,9,17] has achieved fruitful results, such as the reference parameters of the Brazilian disc failure test and PFC2d simulation [9], the laws of various mechanical parameters, the coupling relationship with fluids [12,18], the relationship between the particle size and equilibrium position for the particles in the fluid channel [19], the deformation characteristics of geological bodies, and disaster simulations [20], there are many mechanical interactions between cracks and the particles [5,6]. The method measures particle velocity in water and multiphase flow by particle image velocimetry [21,22]. These studies simulate the soil's circular particles and apply boundary conditions to the circular soil to obtain physical and mechanical results under deformation and stress conditions. The combination of multiple spherical particles can simulate natural soil particles. For example, Sun [12] successfully used different quantities and radii particles to simulate various gravel. Figure 1 shows the connection and composition simulation process from single/circular/ball particles to soil/rock. Many single particle/circular particle/ball particles (Figure 1a) made up a joint particle/crystal particle (Figure 1b). The connection method between single particles is bonded and stacked ( Figure 1c). Moreover, many joint particles ( Figure 1b) make up soils/rock (Figure 1d). The joint particles (Figure 1b) are polygonal and not circle/ball/sphere. DEM uses particles to simulate soils. The simulation is successful when some parameters and results are consistent or the same. In soil mechanics, two crucial indicators of soil mass are particle size and particle gradation. Past research [23] does not concern the particle gradation, or it directly takes the circle diameters as the particle size [24], or the particle radius is adjusted slightly with time stepping [25]. Thus, the past research chose the single particle's ball diameter ( Figure 1a) and did not choose one of the joint particles (Figure 1b) to calculate the particle gradation.
This past simulation has two problems: First, a circular particle (Figure 1a) is not equivalent to a soil particle (Figure 1b). The joint particles are better than the circular particles in representing the soils/rock (Figure 1d). Second, a circular particle's particle size ( Figure 1a) and that of a joint particle (Figure 1b) are different. The particle diameter ( Figure 1a) cannot be used for gradation calculation because actual soil particles are joint particles ( Figure 1b); they are not a circle. The joint particles, not a single circular particle, go through the sieve hole during the soil particles' screening experiment. The joint particle size, not the circular particle size, should be chosen to calculate the gradation. The joint particle (Figure 1b) is the actual particle of soils/rock. This paper aims to solve the above problems, proposes a discrete element, i.e., a joint particle model for soil simulation, and develops a joint particle size calculation method. We used this model and calculation method to study the soil porosity at some depths. Obtain the characteristics of the elastic modulus, porosity, and Poisson's ratio with the depth. The circular discrete element is presented in this research and designed to simulate a soil material more realistically.

Joint Model for Soils
We used many joint soil particles to simulate soil. The paper aims to realize that joint soil particles are cemented and bonded by circular particles to form soil. The number of circles creating a joint soil particle is not the focus of this paper, although scientists need further analysis at this scale. The number was 3-6 in this research. These ball particles connect. One ball particle needs to connect at least one of the other particles. If a ball particle does not relate to any other ball particles, it is dangling and hanging, and it does not belong to this joint particle.
In the numerical simulation test of soil mechanics, we restricted the two sides of the rectangular soil by a wall, and applied pressure through the wall in the soil's upper part. The mechanical interaction between the soil particles ( Figure 2a) and between the particles and the wall plane (Figure 2b,c) involves two effects: pressing and shear sliding translation, simplified into the mechanical model shown in Figure 2. Figure 2d offers the connection and contact between balls connected with two joint particles.
Pressing between two single ball particles means that the distance between the two balls is less than the sum of their radii (Figure 2a). When the distance is greater than or equal to the two radii's sum, the two balls' mechanical action disappears. The interaction between a single circular particle and the wall plane is similar (Figure 2b). The mechanical relationship is F n = K n U n , F s = K s U s , where F n and F s are the normal force and tangential force between two particles or between a particle and a wall plane, respectively. U n is the compression displacement between particles or between a particle and the wall plane. K n is the elastic stiffness. U s is the sliding shear displacement between particles or a single particle and a wall plane. K s is the shear stiffness of this shear.
O j1 , O j2 , O j3, and O j4 are the 1st, 2nd, 3rd, and 4th circular particles in the jth joint particle (Figure 2d). O i1 is the first circular particle in the ith joint particle. The relationship between the distance of two circles' centers and the sum of the two radii is needed to judge whether two joint particles are in contact; to calculate the coincidence degree of mutual extrusion, R im is the mth circular particle radius composed of the ith joint particle (as shown in Figure 2d). R jn is the radius of the nth circular particle composed of the jth joint particle. D is U in Figure 2a. When D ijmn > R im + R jn , two joint particles do not contact each other. When D ijmn = R im + R jn , two joint particles just contact. When D ijmn < R im + R jn , two joint particles squeeze against each other. m = 1, 2, . . . , n = 1, 2, . . . , the maximum value of m and n is the number of circular particles in the joint particle of the ith and jth joint particle. force between two particles or between a particle and a wall plane, respectively. U n is the compression displacement between particles or between a particle and the wall plane. K n is the elastic stiffness. U s is the sliding shear displacement between particles or a single particle and a wall plane. K s is the shear stiffness of this shear.

Figure 2.
Mechanical relationship between ball particles and walls (a). mechanical model of interaction between single ball particles; (b) mechanical model of interaction between particle and wall plane. (c) the relationships. R1, R2 are the radii of single particles; K n is the elastic stiffness, and K s is the shear stiffness. The black shaded area overlaps the ball particles between the ball particles and the wall plane. Three joint particles are red, green, and black, which contact each other. There is a wall plane at the bottom. 1-a, 2-a and 3-a are three contact points between two particles: the contact relationship in (a). 4-b is the contact point between particle and wall, the contact relationship in (b); (d) determination of distance, contact, and extrusion of joint particles.
Oj1, Oj2, Oj3, and Oj4 are the 1st, 2nd, 3rd, and 4th circular particles in the jth joint particle (Figure 2d). Oi1 is the first circular particle in the ith joint particle. The relationship between the distance of two circles' centers and the sum of the two radii is needed to judge whether two joint particles are in contact; to calculate the coincidence degree of mutual extrusion, Rim is the mth circular particle radius composed of the ith joint particle (as shown in Figure 2d). Rjn is the radius of the nth circular particle composed of the jth joint particle. D is U in Figure 2a. When Dijmn > Rim + Rjn, two joint particles do not contact each other. When Dijmn = Rim + Rjn, two joint particles just contact. When Dijmn < Rim + Rjn, two joint particles squeeze against each other. m = 1, 2, ..., n = 1, 2, ..., the maximum value of m and n is the number of circular particles in the joint particle of the ith and jth joint particle. Mechanical relationship between ball particles and walls (a). mechanical model of interaction between single ball particles; (b) mechanical model of interaction between particle and wall plane. (c) the relationships. R 1 , R 2 are the radii of single particles; K n is the elastic stiffness, and K s is the shear stiffness. The black shaded area overlaps the ball particles between the ball particles and the wall plane. Three joint particles are red, green, and black, which contact each other. There is a wall plane at the bottom. 1-a, 2-a and 3-a are three contact points between two particles: the contact relationship in (a). 4-b is the contact point between particle and wall, the contact relationship in (b); (d) determination of distance, contact, and extrusion of joint particles.

Joint Particle Size: Rotation Calculation Model
Some circular balls ( Figure 3a) are connected to form a joint soil particle ( Figure 3b). Some joint particles cluster into soils, which have a specific porosity. Therefore, the circular ball is the elementary unit. These balls, making up a joint particle, have a high connection strength with the other balls in the same soil particle. The stability between balls is high enough to prevent their connections from being broken when soil is under stress to a great extent. More types of joint soil particles are in Figure 3c. Some circular balls (Figure 3a) are connected to form a joint soil particle (Figure 3b). Some joint particles cluster into soils, which have a specific porosity. Therefore, the circular ball is the elementary unit. These balls, making up a joint particle, have a high connection strength with the other balls in the same soil particle. The stability between balls is high enough to prevent their connections from being broken when soil is under stress to a great extent. More types of joint soil particles are in Figure 3c. (b) the maximum width rectangle covering the joint particle; (c) some joint particles with its own the maximum width rectangle; (d) joint particles passing through the sieve by rotation at the maximum width of its rectangle; the solid ball is the elementary ball; the solid black block is the soil particle. Soils consist of many joint soil particles; joint soil particles consist of solid balls. Width and height denote the joint soil particle's width and height, respectively, where Height  Width.
The soil comprises many joint soil particles (Figure 1b), mainly measured and evaluated by the gradation curve. The cumulative percent by the weight of soil passing through a given sieve is the percent of finer. The literature [26,27] showed the gradation curve's details. Choose the uniformity coefficient and the gradation coefficient according to the gradation curve. The uniformity coefficient is Cu = D60/D10; the gradation coefficient is Cc = (D30) 2 /(D60D10), where D60, D30, and D10 are the diameters through which 60%, 30%, and 10% of the total soil mass pass. One must obtain the particle-size distribution of the simulated soil particles constituting the soil to determine the parameters of D60, D30, D10, Cu, and Cc.
The joint particle size is not the size of the balls (the ball diameter); instead, the particle size is the minimum sieve size that passes through. Furthermore, the soil particles' Figure 3. Joint particle size and screening process. (a) Four ball particles composite a joint particle; (b) the maximum width rectangle covering the joint particle; (c) some joint particles with its own the maximum width rectangle; (d) joint particles passing through the sieve by rotation at the maximum width of its rectangle; the solid ball is the elementary ball; the solid black block is the soil particle. Soils consist of many joint soil particles; joint soil particles consist of solid balls. Width and height denote the joint soil particle's width and height, respectively, where Height ≥ Width.
The soil comprises many joint soil particles (Figure 1b), mainly measured and evaluated by the gradation curve. The cumulative percent by the weight of soil passing through a given sieve is the percent of finer. The literature [26,27] showed the gradation curve's details. Choose the uniformity coefficient and the gradation coefficient according to the gradation curve. The uniformity coefficient is C u = D 60 /D 10 ; the gradation coefficient is C c = (D 30 ) 2 /(D 60 D 10 ), where D 60 , D 30, and D 10 are the diameters through which 60%, 30%, and 10% of the total soil mass pass. One must obtain the particle-size distribution of the simulated soil particles constituting the soil to determine the parameters of D 60 , D 30 , D 10 , C u and C c .
The joint particle size is not the size of the balls (the ball diameter); instead, the particle size is the minimum sieve size that passes through. Furthermore, the soil particles' particle size distribution is determined by a sieving test, in which soil particles are passed through a standard sieve, as shown in Figure 3d. The particles passing through the mesh are smaller than the sieve's pore diameter. As a result, the particles remaining on the sieve have larger diameters than the pore diameter of the sieve.
In calculating the particle sizes of joint particles, we considered many rectangles to cover the joint particles. These rectangles cover an entire joint particle and were tangent to the particle boundary. The rectangle with the smallest width (rectangle in Figure 3a-d) among all the rectangles covering the particles was determined. Next, we sieved the joint particles according to Figure 3d; the particles could pass through the sieve only vertically according to the minimum width. If the minimum width of the rectangle was precisely equal to the sieve's diameter, the soil particles could pass through the sieve, and the minimum width is the particle size of the joint soil.
This research proposes a rotation calculation model, as shown in Figure 4. We used the model in two types of conditions. The first condition involved joint soil particles composed of circles, as shown in Figure 4a. The points on the ball boundary had explicit numerical coordinates. Therefore, we only chose boundary points for calculation. The second condition concerned the pixel joint soil particle, as shown in Figure 4b. The dots on the soil particles are pixels, and the joint particle is the picture with pixels. In this method, we chose boundary points and we chose inside points for calculation.  The programming flowchart is shown in Figure 4d. The result w1/w2 and β. Here, w1 and w2 are the minimum width of the joint soil particles and are less than or equal to the sieve diameter in Figure 3d. β = βi is the rotation angle from the initial position when w = w1(βi).
Therefore, the particle size calculation in this paper is a process of discovering the minimum width rectangle. The simulation and indoor screening tests are well unified, establishing a foundation for simulating more realistic soil.

Porosity Estimation of Overlapping Particles: Pixel Counting Method
The porosity of rock and soil in hydrogeology is an important physical parameter that affects and determines the permeability coefficient [28,29]. The porosity should be first calculated to obtain the relationship between soil porosity and pressure. There is no Regarding the joint soil particles composed of circles, the joint particle size is determined by a rotating coordinate system, as shown in Figure 4a. The point coordinates (x, y) on the ball particles were determined, after which we found the new coordinate point (x β , y β ) after the coordinate system rotated an angle β (Figure 4c,d). The relationship between the coordinates of the two coordinate systems is as follows: For a given β in Figure 4c,d, we found the max(x β ), min(x β ), max(y β ) and min(y β ) of the coordinates of all points on the boundaries of all circles, where the two side lengths of rectangles were w 1β = max(x β ) − min(x β ) and w 2β = max(y β ) − min(y β ). When β∈(0, π), min(w 1β ) and min(w 2β ) can be obtained. These two values should be equal, i.e., min(w 1β ) = min(w 2β ) = w. Here, w is the minimum width of the joint soil particles and reflects the joint soil particles' particle size. P(x, y) includes many points. P(x, y) is any point on the joint particles' boundary as shown in Figure 4c.
Similarly, if there are blue particles with arbitrary shapes, as shown in Figure 4b, they comprise blue pixels and other white pixels set to either blue RGB (0, 0, 255) or white RGB (255, 255, 255). We produced any color by blending and adding red, green, and blue (RGB). For example, RGB (200, 225, 255) means that the red color is 200, the green color is 225, and the blue color is 255. The blue color of the two colors, RGB (0, 0, 255) or white RGB (255, 255, 255), has the same number, 255. The first number of RGB, the red color position, was chosen for calculation. Of course, we chose the second number of RGB. The first RGB of the blue pixel position is 0, while the first RGB of the white pixel position is 255. We only changed the 255 into 1. Finally, we made Figure 4b. All points have a value of 1 inside the soils, and each point has corresponding coordinates. According to the calculation method shown in Figure 4b, after rotating the coordinate system by an angle β, calculate the widths w 1β and w 2β corresponding to the minimum rectangle covering a soil particle. When β∈(0, π), min(w 1β ) and min(w 2β ) can be obtained, the values of which should be equal, i.e., min(w 1β ) = min(w 2β ) = w, which is the particle size of the pixel soil particle. According to the image scale, obtain the joint soil particles' actual particle sizes.
The programming flowchart is shown in Figure 4d. The result w 1 /w 2 and β. Here, w 1 and w 2 are the minimum width of the joint soil particles and are less than or equal to the sieve diameter in Figure 3d. β = β i is the rotation angle from the initial position when w = w 1(βi) .
Therefore, the particle size calculation in this paper is a process of discovering the minimum width rectangle. The simulation and indoor screening tests are well unified, establishing a foundation for simulating more realistic soil.

Porosity Estimation of Overlapping Particles: Pixel Counting Method
The porosity of rock and soil in hydrogeology is an important physical parameter that affects and determines the permeability coefficient [28,29]. The porosity should be first calculated to obtain the relationship between soil porosity and pressure. There is no good way to get the porosity of discrete particles in past research [25,30,31]; more researchers choose to avoid or ignore the issue and to not determine the porosity. There was some overlapping areas (the shaded region shown in Figure 2b; B1, B2, B3, and C1 shown in Figure 5a). The porosity is the projected area's rate, not the sum area of all balls in the example area. The porosity is n and 1-ψ, not 1-Φ, shown in Figure 5a. The fundamental difficulty is that the overlapping areas are hard to evaluate and calculate.
This research proposes a more accurate method (as shown in Figure 5) for calculating the soil particle area/volume to solve this problem: The pixel counting method. The soil is formed by accumulating soil particles by gravity (Figure 5b). The soil particles are colored (Figure 5c This research proposes a more accurate method (as shown in Figure 5) for calculating the soil particle area/volume to solve this problem: The pixel counting method. The soil is formed by accumulating soil particles by gravity (Figure 5b). The soil particles are colored (Figure 5c

Elastic Modulus and Poisson's Ratio
We applied pressure to the soils in two directions at specific burial depths. According to the reference [12], we obtained the elastic modulus and Poisson's ratio. We applied Fx and Fy forces in both directions; the stress is given by σ x and σ y , respectively. Fx and Fy are the force applied to the sample in x and y directions. σ x and σ y are the stress on the sample in xand y-direction. The strain formulas are ε x = (σ x − υσ y )/E and ε y = (σ y − υσ x )/E. Because the problem is a plane strain problem, ε x = 0 and σ x = υσ y . Then, σ x + ∆σ x = υ(σ y + ∆σ y ), and, thus, ∆σ x = υ∆σ y . ε x and ε y are the strain on the sample in xand y-direction. υ is Poisson's ratio. E is the elastic modulus.
We applied a slight increase in the force ∆σ y to the soil sample when soils reached stability. ∆σ x can be obtained after the soils reach stability again under ∆σ y . In this way, υ can be obtained from ∆σ x = υ∆σ y . Every pair of ∆σ y and ∆σ y values can correspond to a new value ε y + ∆ε y . We obtained ∆ε y from the soil deformation. Thus, we obtained the elastic modulus E according to ∆ε y = (∆σ y − υ∆σ x )/E.

Joint Particle Size
According to the calculation method in the paper, the crystal particle in Figure 1 is simulated in Figure 6. As the number of circular particles increases in the joint particle, the joint particle shape can better fit the shape of the simulated particles. The crystals of 1, 2, 3, 4, and 5 in Figure 1 can be composed of many circular particles; these circular particles overlap and bond with each other to some extent. As shown in Figure 6, we calculated the particle size according to the method in Section 2.2. From the results, the way presented in this paper is very robust.  Choose crystal 1 in Figure 1 to calculate the particle size. Then, rotate the particle by an angle to find the length of the two rectangular sides. The long side is the height, and the short side is the width. According to the rotation angle and the two lengths, draw Figure 7.

Particle Gradation and Porosity under Pressure
The weight of the upper soil compresses the soil at a specific depth. The circular particles in contact with each other once again overlap. Figure 2 shows the shadow region, and the pore sizes become small. We established the compressed result image and obtained the porosity according to the image pixel count method presented in Section 2.3.
Over a long geological history, the soil has gradually accumulated small particles. As a result, the thickness and upper pressure have steadily increased. The upper pressure imposed on a specific soil depth is related to gravity's effect on the upper soil. P = gh, where  = s(1 − n), P is the upper pressure,  is the soil density, s is the soil particle density, g is the gravitational acceleration, h is the depth, and n is the soil porosity. Furthermore, we obtained the relationship between the upper pressure and porosity.
We chose many physical parameters ( Figure 2) [33,34]. The ball density  is 2300 kg/m 3 . The normal contact stiffness K n is 4.4 × 10 7 N/m, and the shear contact stiffness K s is 2.2 × 10 7 N/m. The ratio of the two stiffness K s /K n is 0.5. The gravitational acceleration g is 9.81 m/s 2 , the ball radius r is 2-3 mm, and the friction coefficient Fr is 0.5. This paper presents an example calculation. Because the calculation needs to consume a lot of calculation time, we chose the joint particles composed of 3-5 circular particles for research in this paper. We simulated other, more circular, particles in the same method and steps. The research process was as follows: (1) Generate two thousand joint particles (Figures 1b and 3b,c). The two thousand joint particles comprise circular particles (Figure 1a and balls in Figure 3) sized and located randomly according to a specific distribution method, such as the normal, index, and lognormal, uniform distribution. We generated the circle center by the normal distribution in the following example, and the uniform distribution chose the radii. (2) A certain number of soil particles (Figures 1b and  3b,c) were randomly selected from the 2000 soil particles and were allowed to accumulate In this paper, the particle size is the width of the minimum width rectangle, the minimum value of the solid red curve in Figure 7. The rectangle with the minimum width is (5.441 mm, 3.707 mm). The rotation angles are 0.30 and 1.87, that is, 17.189 • and 107.143 • with an error of 0.046 • , which is acceptable. According to previous studies, the average size of the rectangle length is the dotted line. The minimum value is 4.574 mm. However, the green dotted line and the solid red line are different, and the two curves do not coincide. Therefore, the minimum values are 3.707 mm and 4.574 mm, which are not equal. According to Figure 3 and Section 2.2, the particle size in this paper is scientific and reasonable; the calculation process conformed to the test environment and conditions. Therefore, it is inappropriate to take the average value of the sum of the length of the previous value as the particle size.
The length and width of the rectangle changed simultaneously with the rotation angle, as shown in Figure 7. It shows the rule of periodic change with period π/2. In the calculation, taking the value (0, π/2), not (0, π), to calculate the rotation range can meet the requirements. More joint particle lengths and widths have been calculated and a database built. Refer to the literature [32].

Particle Gradation and Porosity under Pressure
The weight of the upper soil compresses the soil at a specific depth. The circular particles in contact with each other once again overlap. Figure 2 shows the shadow region, and the pore sizes become small. We established the compressed result image and obtained the porosity according to the image pixel count method presented in Section 2.3.
Over a long geological history, the soil has gradually accumulated small particles. As a result, the thickness and upper pressure have steadily increased. The upper pressure imposed on a specific soil depth is related to gravity's effect on the upper soil. P = ρgh, where ρ = ρ s (1 − n), P is the upper pressure, ρ is the soil density, ρ s is the soil particle density, g is the gravitational acceleration, h is the depth, and n is the soil porosity. Furthermore, we obtained the relationship between the upper pressure and porosity.
We chose many physical parameters ( Figure 2) [33,34]. The ball density ρ is 2300 kg/m 3 . The normal contact stiffness K n is 4.4 × 10 7 N/m, and the shear contact stiffness K s is 2.2 × 10 7 N/m. The ratio of the two stiffness K s /K n is 0.5. The gravitational acceleration g is 9.81 m/s 2 , the ball radius r is 2-3 mm, and the friction coefficient Fr is 0.5. This paper presents an example calculation. Because the calculation needs to consume a lot of calculation time, we chose the joint particles composed of 3-5 circular particles for research in this paper. We simulated other, more circular, particles in the same method and steps. The research process was as follows: (1) Generate two thousand joint particles (Figures 1b and 3b,c). The two thousand joint particles comprise circular particles (Figure 1a and balls in Figure 3) sized and located randomly according to a specific distribution method, such as the normal, index, and lognormal, uniform distribution. We generated the circle center by the normal distribution in the following example, and the uniform distribution chose the radii. (2) A certain number of soil particles (Figures 1b and 3b,c) were randomly selected from the 2000 soil particles and were allowed to accumulate under the effect of gravity. (3) We applied the pressure to the upper part of the generated soil. The magnitude of the pressure corresponded to the effect of gravity on the soil in this depth range. (4) We calculated the soil porosity under this pressure. We obtained the relationship between the upper pressure and the porosity. In the range of x ∈ (−0.02 m, 0.1 m), y ∈ (−0.02 m, 1.5 m), 400, 800, 1200, 1600, and 2000 soil particles were selected, and C u and C c were calculated according to Section 2.2. The results are listed in Table 1 and shown in Figure 8. Table 1. C u and C c of soils made up of two types of particles. The particle min-width in Figure 8a is the joint particle (Figures 3 and 4) size of soils, the minimum rectangle's minimum width covering the joint soil particle. The particle min-width in Figure 8b is the single ball's radius (Figures 3 and 4).

Number of Particles
C u and C c in Table 1 of soils were made up of two types of particles. D 60 , D 30, and D 10 were the sieve diameter, with 60%, 30%, and 10% of the total soil mass passing through the sieve (unit: mm). C u = D 60 /D 10 was the uniformity coefficient. C c = (D 30 ) 2 /(D 60 D 10 ) was the gradation coefficient. C cc and C cs were C c of the combined particles and the single ball, respectively. C uc and C us are the C u of the combined particles and the single ball, respectively. | | denotes the absolute value. C u = D 60 /D 10 is the uniformity coefficient. C c = (D 30 ) 2 /(D 60 D 10 ) is the gradation coefficient. C u , C c are important parameters in particle, soil, and sand research, and engineering [35]. We used the two parameters to distinguish between well graded and poorly graded coarse-grained soil using laboratory tests of the grain size distribution [36]. The two parameters satisfy the requirements in engineering projects or science research. We can say that the materials were qualified [37]. Figure 8a is a gradation curve for the joint particles, similar to some of the literature's gradation curves [38,39]. Figure 8a and Table 1 show that the five gradation curves were not different, and that the difference between C u and C c was not significant. Therefore, the gradation curve was not closely related to the number of selected soil particles. Figure 8b is a gradation curve for ball particles and not joint particles. Figure 8b and Table 1 show that the five gradation curves were almost coincident. Figure 8a,b shows that the particle size in Figure 8a was larger than the particle size in Figure 8b and that the corresponding gradation curves were significantly different. Therefore, we did not take the ball particle gradation curve as the joint soil particles' gradation curve, which is different from previous research methods [12]. boundary pressure and porosity. Next, we obtained the elastic modulus and Poisson's ratio according to the calculation procedure mentioned in Section 2.4, as shown in Figure  11. We fit the following formula to the porosity and depth in Figure 10.

Depth n a b c =  +
(2) a = 0.311, b = 0.973, c = 0.0378; the correlation coefficient R was 0.998. We pressurize the sample from the top and restricted the displacement on both sides. The pressures above and on both sides were σ1 and σ2, respectively. Hσ1 and Hσ2 are the pressure conversion depths of σ1 and σ2, respectively, where Hσ1 = σ1/(g) and Hσ2 = σ2/(g); h is the depth;  is the soil density; g is the acceleration due to gravity. This gradation curve of joint particles in Figure 8a has not been calculated and displayed in previous research. The reason is that the particle size of joint soil particles could not be well estimated. We obtained the particle size of the soil particles via the calculation method of this paper. First, we determined the gradation curve of the granular soil. After compression, the range of images taken was x∈(−0.02 m, 0.1 m), y∈(−0.02 m, 0.1 m), and the depth ranged from 10 m to 90 m, with an interval of 1 m. We recorded the force on the upper, lower, left, and proper boundaries and the sample's height with the constant increase in the upper pressure. The compressed soil picture calculated the porosity using the pixel recognition technology discussed in Section 2.3. Figures 9 and 10 show the boundary pressure and porosity. Next, we obtained the elastic modulus and Poisson's ratio according to the calculation procedure mentioned in Section 2.4, as shown in Figure 11.
We fit the following formula to the porosity and depth in Figure 10.
a = 0.311, b = 0.973, c = 0.0378; the correlation coefficient R was 0.998. We pressurize the sample from the top and restricted the displacement on both sides. The pressures above and on both sides were σ 1 and σ 2 , respectively. H σ1 and H σ2 are the pressure conversion depths of σ 1 and σ 2 , respectively, where H σ1 = σ 1 /(ρg) and H σ2 = σ 2 /(ρg); h is the depth; ρ is the soil density; g is the acceleration due to gravity. Figure 9 shows that the change law of the upper and lower pressures had a good consistency. H σ1 was more significant than H σ2 , and both forces were greater than the natural depth h. The slope of Line h is 1. The slopes of Line H σ2 and Line H σ1 gradually increase with the depth increasing. After 70 m, the slopes of Line H σ2 and Line H σ1 tend to decrease, at approximately 80 m, Line H σ2 begins to fall below the h line.   The curve n in Figure 10 shows that the porosity gradually decreased as the depth increased. The porosity gradually reduced with the depth increasing from the image analysis of the soil in P1, P2, P3, P4, and P5. Figure 10 shows the following: (1) when the depth was small, some interconnected pores existed between the soil particles, and the porosity was large; (2) when the depth was significant, the particles were superimposed on each the porosity curve tended to be stable; from about 80 m, the porosity tendd to increase slightly. The soil particles' elastic modulus and Poisson's ratio gradually increased as the depth increased ( Figure 11). The elastic modulus and Poisson's ratio underwent some upward and downward fluctuations. At a depth of about 48 m, the elastic modulus had more fluctuations.
In contrast, Poisson's ratio maintained a relatively steady increase. The fluctuation of the elastic modulus was more significant than that of Poisson's ratio. With the depth increase, Poisson's ratio fluctuated when about 62 m.

Discussion
As the soil depth increased, the upper pressure increased, the distance between the soil particles gradually decreased, many soil particles appear to overlap (from the progression order in P1, P2, P3, P4, and P5 in Figure 10), and the pore sizes between the soil particles decreased. Therefore, as the depth of burial increased, the porosity (shown in Figure 10) gradually decreased.
With more overlap between soil particles, the soil's ability to resist external force deformation was enhanced. Under the same external pressure and stress σ, the deformation and strain ε became smaller. According to the formula E = σ/ε, the elastic modulus ( Figure  11) gradually increased as the depth increased. The parameters in Equation (2) can preliminarily clarify the physical meaning: Parameter c is the porosity when the Depth → ∞, which indicates that no pores exist between the particles, the porosity only pertains to the pores inside the material, and this part is related to the microscopic morphology inside the material. When depth → 0, a + c is the porosity; the soil porosity is n0 in an entirely pressure-free and loose state. At that time, a = n0 − c is related to the loose state of the material. B represents the degree of porosity increase with the depth increasing. B has a direct relationship with the elastic modulus E of the material and the density  of the particles and is, thus, related to the material's mechanical parameters.
At 70 m, the porosity tended to be stable. From 80 m, the porosity curve increased slightly. The fluctuation occurred because the particles detached from the soil [10]. In calculating the PDEM, we applied the external force to the soil particles by moving the wall toward the soil [44]. As the wall moved toward the soil, it compressed the soil. As a result, the interaction force between the soil particles and the wall gradually increases. According to Figure 2b and the formula F n = K n U n , U n must be less than R. As the pressure increased, Curve n and curve ψ in Figure 8 are results obtained by this research, while curve Φ is the result ignored in past research. The directional area accumulation of balls is not the projected area; the porosity n was 1-ψ, not 1-Φ from Figures 5a and 10.
The curve n in Figure 10 shows that the porosity gradually decreased as the depth increased. The porosity gradually reduced with the depth increasing from the image analysis of the soil in P1, P2, P3, P4, and P5. Figure 10 shows the following: (1) when the depth was small, some interconnected pores existed between the soil particles, and the porosity was large; (2) when the depth was significant, the particles were superimposed on each other, the connected pores decreased in number or even disappeared, the pore sizes gradually decreased, and the porosity decreased. This characteristic is consistent with the literature's conclusion that the density increases [40][41][42][43]. Thus, the density of soil after compression increased, while the porosity decreased. From a depth of approximately 72 m, the porosity curve tended to be stable; from about 80 m, the porosity tendd to increase slightly. The soil particles' elastic modulus and Poisson's ratio gradually increased as the depth increased ( Figure 11). The elastic modulus and Poisson's ratio underwent some upward and downward fluctuations. At a depth of about 48 m, the elastic modulus had more fluctuations.
In contrast, Poisson's ratio maintained a relatively steady increase. The fluctuation of the elastic modulus was more significant than that of Poisson's ratio. With the depth increase, Poisson's ratio fluctuated when about 62 m.

Discussion
As the soil depth increased, the upper pressure increased, the distance between the soil particles gradually decreased, many soil particles appear to overlap (from the progression order in P1, P2, P3, P4, and P5 in Figure 10), and the pore sizes between the soil particles decreased. Therefore, as the depth of burial increased, the porosity (shown in Figure 10) gradually decreased.
With more overlap between soil particles, the soil's ability to resist external force deformation was enhanced. Under the same external pressure and stress σ, the deformation and strain ε became smaller. According to the formula E = σ/ε, the elastic modulus ( Figure 11) gradually increased as the depth increased. The parameters in Equation (2) can preliminarily clarify the physical meaning: Parameter c is the porosity when the Depth → ∞, which indicates that no pores exist between the particles, the porosity only pertains to the pores inside the material, and this part is related to the microscopic morphology inside the material. When depth → 0, a + c is the porosity; the soil porosity is n 0 in an entirely pressure-free and loose state. At that time, a = n 0 − c is related to the loose state of the material. B represents the degree of porosity increase with the depth increasing. B has a direct relationship with the elastic modulus E of the material and the density ρ of the particles and is, thus, related to the material's mechanical parameters.
At 70 m, the porosity tended to be stable. From 80 m, the porosity curve increased slightly. The fluctuation occurred because the particles detached from the soil [10]. In calculating the PDEM, we applied the external force to the soil particles by moving the wall toward the soil [44]. As the wall moved toward the soil, it compressed the soil. As a result, the interaction force between the soil particles and the wall gradually increases. According to Figure 2b and the formula F n = K n U n , U n must be less than R. As the pressure increased, K n remained stable during the calculation; U n gradually increased. When U n > R, the black shaded area in Figure 2b exceeded half the circle; the circle center passed through the wall. The force between the ball and the wall changed to F n = K n (2R − U n ), and the direction changed by 180 • . The ball immediately passed through and away from the wall.
The elastic modulus and Poisson's ratio volatility in Figure 11 are greater than the porosity volatility in Figure 10. Under the above pressure, the particles are superimposed and elastically deformed. When the particles escape, the force on the boundary changes immediately. According to the formula ∆σ x = υ∆σ y and ∆ε y = (∆σ y − υ∆σ x )/E [12], the elastic modulus and Poisson's ratio change drastically. However, space opened up by particles' ejection is immediately occupied by the elastic release and deformation recovery of other elastically deformed particles. As a result, the superimposed elastic area between the particles (as shown in Figure 2a,b) was diminished. The release of the particle superposed area makes up for the space left by escaping particles. The particle area in the soil does not change much; thus, the soil porosity does not change significantly.
The force in the y direction produced a strain ε y in y direction. Similarly, the force in the x direction produced a strain ε x in the y direction. The y direction was the primary stress direction. Poisson's ratio is the ratio of the lateral strain ε x to the active strain ε y , i.e., υ = ε x /ε y . When the pressure in the active direction was slight, the specimen was not prone to the lateral strain; Poisson's ratio was relatively small. However, as the pressure increased and the particles were pressed and superimposed, the entire soil became more compact. After the strain occurred in the active direction, the lateral direction was more prone to strain. As a result, Poisson's ratio increased. Therefore, Poisson's ratio increased with the burial depth and the upper load applied to the soils (as shown in Figure 11).

Conclusions
This paper was founded on the model and method, the rotation calculation model for the particle size and pixel counting method for the porosity. Further, we obtained the particle gradation of joint particles. The two models overcame past research shortcomings, e.g., considering soil particles as only circles and balls [45][46][47][48], ignoring the particle size and particle gradation of joint soil particles and the actual porosity. This research enables the development of more realistic discrete elements and helps to simulate the more complex rock and soil materials. In particular, this study provides a good reference for the transport of foundation sediment in porous media [49] and containing faults [29]. The founded calculation model combined with the discrete element force calculation. The relationship between porosity, elastic modulus, Poisson's ratio, and soil pressure (the buried depth) was studied. Soil, composed of soil particles, was compressed by an external force. The compactness, Poisson's ratio, and the elastic modulus increased, while the porosity decreased, consistent with previous research results [50,51].
When soil particles were ejected from the soil, and the depth was more than the sensitive depth values, Poisson's ratio, the elastic modulus, and the porosity had significant fluctuations, but the elastic modulus changed significantly. Therefore, the elastic modulus is a good index and is the physical parameter for primary monitoring in the integrity monitoring of underground petroleum and groundwater reservoirs.
This research focused on the model, calculating the particle size of soil particles and studying the relationship between compression porosity, the elastic modulus, Poisson's ratio, and depth. We did not consider the ball radius size and the number when simulating soil particles, and, therefore, did not consider the scale of the soil particle size. Simulating soil particles using different radii and different numbers of balls may yield different final results and be considered in subsequent studies.